function twoshiftpdepsmiddle(N,D1,D2,M,eps1,eps2)
% twoshiftpdepsmiddle(N,D1,D2,M,eps1,eps2)
% 	Parameters:
% 		N:	number of data points to include in the time series
% 		D1:	first number of seconds to displace time series over
% 		D2:	second number of seconds to displace time series over
% 		M:	kernel size
% 		eps1:	initial amount to perturb by
% 		eps2:	second amount to perturb by
% 	Outputs:
% 		Only the plots this code produces.
% 	Purpose of this code:  This code is a modification of the code twoshiftpdeps.  Here the interpolation is centered at the mean value of the 
% 		time series, as opposed to the inital value.
hold off;

K = 1:N;
ts = whitenoise(N,'w',1);
ts = lowpass(ts);
%ts = testsine(K);
eps = eps1 + eps2;
D = D1+D2;



D1M = D1+((N./2)).*eps1;
D2M = D2+((N./2)).*eps2;
DM = D+((N./2)).*eps;

[psorig,freq] = mypsd(ts);
%% shift timeseries by S1 and then by S2
F11 = FDtesteps(ts,D1M,M,@sinc);
F12 = FDtesteps(ts,D1M,M,@sincdir);
F13 = FDtesteps(ts,D1M,M,@sinc2dir);


for j = 1:N
    %delta(j) = (D1 + (j - 1).*eps1) - D1;
    delta(j) = (j-1-(N./2)).*eps1;
end

for j = 1:N
    tsint10(j) = F11(j);
    tsint11(j) = F11(j) + delta(j).*F12(j);
    tsint12(j) = F11(j) + delta(j).*F12(j) + (1./2).*(delta(j)^2).*F13(j);
end

for j = 1:N
    Dadj = D1M - (M-1)/2; 
    Dint = round(Dadj); 
    k = (-(M-1)/2):1:((M-1)/2); 
    s = (Dadj-Dint)-k; 
    w = getWindow('blackman3',M);
    p = sinc(s).*w;
    q = sincdir(s).*w;
    q2 = sinc2dir(s).*w;
    
    intf0 = (p);% + eps.*q + (1./2).*(eps).^2.*q2);
    intf1 = (p + delta(j).*q);% + (1./2).*(eps).^2.*q2);
    intf2 = (p + delta(j).*q + (1./2).*(delta(j)).^2.*q2);
    %intf = (p + q);
    n0(j) = sum(intf0);    
    n1(j) = sum(intf1);    
    n2(j) = sum(intf2);
end

tsint10 = tsint10./n0;
tsint11 = tsint11./n1;
tsint12 = tsint12./n2;

tsd1 = tsint11;

F21 = FDtesteps(tsd1,D2M,M,@sinc);
F22 = FDtesteps(tsd1,D2M,M,@sincdir);
F23 = FDtesteps(tsd1,D2M,M,@sinc2dir);

for j = 1:N
    %delta(j) = (D2 + (j - 1).*eps2) - D2;
    delta(j) = (j-1-(N./2)).*eps2;
end

for j = 1:N
    tsint20(j) = F21(j);
    tsint21(j) = F21(j) + delta(j).*F22(j);
    tsint22(j) = F21(j) + delta(j).*F22(j) + (1./2).*(delta(j)^2).*F23(j);
end

for j = 1:N
    Dadj = D2M - (M-1)/2; 
    Dint = round(Dadj); 
    k = (-(M-1)/2):1:((M-1)/2); 
    s = (Dadj-Dint)-k; 
    w = getWindow('blackman3',M);
    p = sinc(s).*w;
    q = sincdir(s).*w;
    q2 = sinc2dir(s).*w;
    
    intf0 = (p);% + eps.*q + (1./2).*(eps).^2.*q2);
    intf1 = (p + delta(j).*q);% + (1./2).*(eps).^2.*q2);
    intf2 = (p + delta(j).*q + (1./2).*(delta(j)).^2.*q2);
    %intf = (p + q);
    n0(j) = sum(intf0);    
    n1(j) = sum(intf1);    
    n2(j) = sum(intf2);
end

tsint20 = tsint20./n0;
tsint21 = tsint21./n1;
tsint22 = tsint22./n2;

tsd2 = tsint21;

ps = mypsd(tsd2);


F1 = FDtesteps(ts,DM,M,@sinc);
F2 = FDtesteps(ts,DM,M,@sincdir);
F3 = FDtesteps(ts,DM,M,@sinc2dir);

for j = 1:N
    %delta(j) = (D + (j - 1).*eps) - D;
    delta(j) = (j-1-(N./2)).*eps;
end

for j = 1:N
    tsint0(j) = F1(j);
    tsint1(j) = F1(j) + delta(j).*F2(j);
    tsint2(j) = F1(j) + delta(j).*F2(j) + (1./2).*(delta(j)^2).*F3(j);
    
end
for j = 1:N
    Dadj = D - (M-1)/2; 
    Dint = round(Dadj); 
    k = (-(M-1)/2):1:((M-1)/2); 
    s = (Dadj-Dint)-k; 
    w = getWindow('blackman3',M);
    p = sinc(s).*w;
    q = sincdir(s).*w;
    q2 = sinc2dir(s).*w;
    
    intf0 = (p);% + eps.*q + (1./2).*(eps).^2.*q2);
    intf1 = (p + delta(j).*q);% + (1./2).*(eps).^2.*q2);
    intf2 = (p + delta(j).*q + (1./2).*(delta(j)).^2.*q2);
    %intf = (p + q);
    n0(j) = sum(intf0);    
    n1(j) = sum(intf1);    
    n2(j) = sum(intf2);
end
tsint0 = tsint0./n0;
tsint1 = tsint1./n1;
tsint2 = tsint2./n2;

tsonce = tsint1;
psa = mypsd(tsonce);

for j = 1:N
    tsy1a(j) = testsine(j-(D1+(j-1).*eps1));
end

E1 = tsd1 - tsy1a;

F31 = FDtesteps(tsy1a,D2M,M,@sinc);
F32 = FDtesteps(tsy1a,D2M,M,@sincdir);
F33 = FDtesteps(tsy1a,D2M,M,@sinc2dir);

for j = 1:N
    %delta(j) = (D2 + (j - 1).*eps2) - D2;
    delta(j) = (j-1-(N./2)).*eps2;
end

for j = 1:N
    tsint30(j) = F31(j);
    tsint31(j) = F31(j) + delta(j).*F32(j);
    tsint32(j) = F31(j) + delta(j).*F32(j) + (1./2).*(delta(j)^2).*F33(j);
end

for j = 1:N
    Dadj = D2M - (M-1)/2; 
    Dint = round(Dadj); 
    k = (-(M-1)/2):1:((M-1)/2); 
    s = (Dadj-Dint)-k; 
    w = getWindow('blackman3',M);
    p = sinc(s).*w;
    q = sincdir(s).*w;
    q2 = sinc2dir(s).*w;
    
    intf0 = (p);% + eps.*q + (1./2).*(eps).^2.*q2);
    intf1 = (p + delta(j).*q);% + (1./2).*(eps).^2.*q2);
    intf2 = (p + delta(j).*q + (1./2).*(delta(j)).^2.*q2);
    %intf = (p + q);
    n0(j) = sum(intf0);    
    n1(j) = sum(intf1);    
    n2(j) = sum(intf2);
end

tsint30 = tsint30./n0;
tsint31 = tsint31./n1;
tsint32 = tsint32./n2;

tsdy2 = tsint31;

tsdy2a = testsine(K-(D+(j-1).*(eps1+eps2)));

E2 = tsdy2 - tsdy2a;

E2hat = tsd2 - tsdy2a;

%figure;
%plot(1:N,E1,'r');

%figure;
%plot(1:N,E2,'b');

%figure;
%plot(1:N,E2hat,'g');




%% First Plot
%loglog(1:length(psorig),psorig,'b');
%hold on;
%err = abs(ps - psa)./abs(psa);
%xnn = abs(ps - psa);
%supp = sqrt(psa./xnn);
%loglog(1:length(ps),err,'g');
%loglog(1:length(ps),supp,'r');
%plot(1:length(tsd1),tsd1,'r');
err2 = abs(tsd2-tsonce);
err3 = err2./abs(tsonce);
%figure;
%semilogy(1:N,err3,'g');
xlim([300,N]);
errps = mypsd(err2);
%plot(1:length(tsd2),err2,'r');
%xlim([300,length(tsd2)]);
%figure;
%loglog(freq,errps,'g');
suppr = sqrt(psorig./errps);
figure;
loglog(freq,suppr);
%hold on;

